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Abstract 

For  purposes  relating  to  the  U.S.  Army’s  need  for  materials  modeling  and 
force  protection,  this  work  provides  justification  for  assigning  effective 
equivalence  between  two  commonly  used  fluid  simulation  methods— 
namely  the  Navier-Stokes  (NS)  and  Lattice  Boltzmann  methods.  The 
Lattice  Boltzmann  Method  (LBM)  has  become  increasingly  popular  as  an 
alternative  approach  to  traditional  NS-based  techniques  for  modeling 
various  incompressible  fluid  flow  applications.  The  LBM  has  recently 
increased  its  range  of  applicability  to  include  numerous  fields  of  interest 
including  those  involving  multiphase  and  thermo-fluid  structure 
interactions.  This  report  documents  a  comparison/validation  effort 
accompanying  the  development  of  a  standard  Lattice  Boltzmann  solver 
with  immersible  moving  boundaries.  The  primary  goal  is  to  validate  the 
model  by  comparing  it  with  various  laminar,  incompressible  flow  cases 
simulated  using  a  finite  volume-based  NS  solver.  Simulations  involving 
four  standard  benchmark  studies  were  analyzed:  (i)  the  flow  through  a 
rectangular  channel,  (2)  the  flow  through  a  lid-driven  cavity,  (3)  the  flow 
over  a  back-step,  and  (4)  the  flow  over  a  stationary  circular  cylinder.  For 
the  specific  applications  and  Reynolds  numbers  simulated,  the  results 
showed  excellent  agreement  between  the  two  cases.  Disparities  were 
observed  only  when  the  theoretical  constraints  of  the  LBM  were  exceeded. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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1  Introduction 

1.1  Background 

The  Lattice  Boltzmann  method  (LBM)  is  becoming  increasingly  popular  as 
an  alternative  approach  to  traditional  techniques  for  modeling  various  flu¬ 
id  flow  applications.  More  traditional  approaches  such  as  the  Navier- 
Stokes  (NS)-based  finite  difference  method  (FDM),  the  finite  element 
method  (FEM),  or  the  finite  volume  method  (FVM)  are  founded  on  the 
discretization  of  macroscopic  equations  at  the  continuum  level.  By  con¬ 
trast,  the  LBM  is  based  on  distribution  functions  describing  the  kinetic  be¬ 
havior  of  particles  that  are  characteristic  of  microscale  and  mesoscale 
simulations  (Griebel;  et  al.  1997;  Al-Jahmany  2004;  Zienkiewicz  et  al. 
2005;  Canuto  et  al.  1988). 

The  use  of  the  LBM  approach  for  numerically  simulating  secondary 
(amorphous)  phase  constituents  (most  often  represented  as  a  laminar,  in¬ 
compressible  fluid)  within  the  context  of  material  science-based  applica¬ 
tions  is,  for  our  purposes,  strategically  based  on  two  primary  motivations. 
The  first  of  these  motivations,  and  possibly  the  most  dominant,  relates  to 
the  omnipresent  problem  of  spatiotemporal  scale  compatibility  and  recon¬ 
ciliation,  wherein  one  attempts  to  make  macroscale  material  decisions 
based  on  microscale  or  nanoscale  performance  criteria.  The  second  moti¬ 
vation  relates  to  concerns  over  mesh  quality  and  the  ability  of  the  method 
being  used  to  accurately  resolve  the  flow  over  complex  shapes. 

With  respect  to  scale  compatibility,  the  LBM  is  derived  from  a  strictly  at¬ 
omistic  (though  probabilistic)  parent  method.  Known  as  the  Boltzmann 
equation,  the  parent  method  describes  the  phenomenological  properties  of 
microscale,  particulate  systems  via  the  use  of  a  phase  space  distribution 
function.  Particularly  attractive,  is  the  fact  that  the  method  allows  for  the 
direct  computation  of  macroscopic  variables  (e.g.,  velocity,  density,  pres¬ 
sure)  through  a  system  of  moment-based  equations,  as  well  as  information 
at  the  microscale,  vis-a-vis  direct  interpretation  of  the  distribution  func¬ 
tion.  The  LBM  derivation  parallels  its  parent  method  in  every  particular, 
except  that  it  allows  for  a  more  simplified  estimation  of  the  probability  dis¬ 
tribution  function  via  a  set  of  discrete  velocities  situated  on  a  structured 
lattice  framework.  Multiscale  applicability  is  thus  also  an  inherent  feature 
within  this  method,  as  it  has  successfully  demonstrated  its  capacity  to 
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solve  a  number  of  problems  at  both  the  micro  and  continuum  scales 
(Raabe  2004;  Sued  et  al.  1993).  The  more  traditional  methods  (e.g.,  Euler, 
NS,  Burnett),  though  arguably  more  ubiquitous  in  their  usage  (owing  to 
their  extended  level  of  maturity,  and  familiarity),  are  nevertheless  derived 
strictly  from  a  continuum-based  assumption.  Apart  from  a  few  exceptions 
(e.g.,  the  use  of  slip  wall  boundary  conditions  to  extend  the  NS  method  in¬ 
to  regimes  of  larger  Knudsen  number  [Kn]),  these  methods  are  quite  lim¬ 
ited  for  purposes  of  spatiotemporal  scaling. 

The  second  argument  in  support  of  using  the  LBM  approach  concerns  the 
issue  of  accurately  resolving  the  flow  field  dynamics  surrounding  complex 
shapes.  Grain  structures  (static,  dynamic,  isotropic,  anisotropic),  structur¬ 
al  asperities,  irregular  flow  channels,  inclusions,  etc.  are  common  features 
found  in  many  material  science  applications.  Most  traditional  methods  re¬ 
ly  on  highly  resolved  (unstructured)  grid  approaches  to  model  these  com¬ 
plex  geometries.  Unfortunately,  these  traditional  methods  are  subject  to 
several  potential  drawbacks,  including: 

•  Computational  costs  due  to  the  large  volume  of  elements  required  to 
adequately  resolve  “small  tight  spaces”  and  areas  of  large  flow  gradi¬ 
ents. 

•  Concerns  over  mesh  quality,  which  inevitably  occur  in  conjunction  with 
the  computation  cost  drawback  just  above  and  include  such  items  as 
negative  element  volumes,  highly  skewed  elements,  and  improper  sur¬ 
face  normal  assignments. 

•  Time  and  resource  costs  due  to  the  excessive  labor  requirements  to  im¬ 
plement  quality  grid(s).  (This  often-underestimated  concern  consumes 
the  vast  majority  of  a  researchers’  time  and  resources,  often  resulting 
in  the  separate  hire  of  a  “meshing  expert”). 

The  LBM  approach  is  effectively  immune  from  these  drawbacks  due  in 
large  part  to  the  development  of  the  immersible  moving  boundary  (IMB) 
method  (Owen  et  al.  2011).  While  the  LBM  does  utilize  a  simple  structured 
regular  lattice,  the  IMB  method  operates  in  a  manner  that  is  completely 
independent  and  distinct  from  the  underlying  mesh  density.  The  IMB 
method  is  based  on  the  adaptation  of  the  LBM  collision  operator  which 
accounts  for  the  limited  area  or  volume  fraction  of  a  lattice  cell,  subject  to 
an  intruding  surface  or  volume.  While  it  is  still  necessary  to  provide  a  reli¬ 
able  method  for  the  computation  of  the  said  area  or  volume  fraction  (de- 
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pending  on  the  prescribed  level  of  accuracy),  there  are  typically  a  variety  of 
proven  methods  that  are  readily  practicable. 

Additional  advantages  of  the  LBM  over  the  more  conventional  approaches 
include  the  following: 

•  The  LBM  is  relatively  simple  to  implement.  Unlike  other  solvers,  LBM 
models  flow  fields  explicitly.  Each  time  step  can  be  split  into  a  stream¬ 
ing  step  and  a  collision  step  involving  no  matrix  inversion  processes. 
The  development  of  the  standard  Bhatnagar-Gross-Krook  (BGK)  colli¬ 
sion  operator  (Higuera  and  Jimenez  1989;  Koelman  1991)  further  eases 
implementation  by  providing  a  dramatically  simplified  alternative  to 
the  cumbersome  integrodifferential  collision  operator  of  the  standard 
Boltzmann  equation. 

•  The  LBM  provides  excellent  scalability  for  large,  parallel  applications. 
All  interactions  in  the  LBM  are  strictly  local.  In  the  collision  step,  the 
updated  distribution  function  at  a  given  lattice  site  is  computed  by  in¬ 
volving  only  the  information  at  the  same  node.  In  the  propagation  step, 
each  lattice  site  must  exchange  information  only  with  its  nearest- 
neighbor  lattice  sites. 

•  The  convection  term  in  the  LBM  is  linear.  In  contrast  to  the  second- 
order  non-linear  convection  term  in  the  NS  equations,  the  linearity  of 
the  LBM  convection  term  results  in  significant  time  savings.  The  LBM 
recovers  the  same  second-order  accuracy  as  the  NS  equations  through 
a  relaxation  process  involving  particle  interactions. 

•  Unlike  NS-based  solutions,  an  explicit  solution  for  the  pressure  dynam¬ 
ics  is  not  required  in  the  LBM  approach.  The  pressure  in  the  LBM  is 
computed  from  the  ideal  gas  law.  The  NS  approach  involves  a  separate 
iterative  solution  involving  the  Poisson  equation.  This  calculation  is 
one  of  the  most  time-consuming  operations  of  the  NS  method. 

As  a  consequence  of  these  recent  advances,  the  LBM  has  extended  its 
range  of  applicability  to  include: 

•  simulations  of  turbulent  flows  (Eggels  1996), 

•  suspension  of  colloidal  particles  (Ladd  2001), 

•  porous  media  (Gunstensen  et  al.  1991), 

•  multiphase  and  multicomponent  flows  (Grunau  et  al  1993;  Shan  and 
Chen  1993),  and 

•  magnetohydrodynamics  (Chen  et  al.  1991). 
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Other  LBM-based  research  efforts  that  have  been  specifically  conducted  at 
ERDC  include: 

•  simulation  of  viscous  flow  through  a  column  of  glass  beads  (Maier  et  al. 

1998), 

•  investigations  involving  the  effects  of  pore  geometry  on  transport 
mechanisms  to  better  characterize  constitutive  relationships  in  multi¬ 
phase  flows  (Maier  et  al.  1996),  and 

•  simulations  involving  the  accuracy  of  the  LBM  as  relating  to  Mach 
number,  spatial  resolution,  boundary  conditions  and  the  lattice  mean 
free  path  (Maier  and  Bernard  1997). 

Additionally,  or  in  many  cases  complementary  to  the  aforementioned  list¬ 
ing,  the  LBM  has  recently  been  developed  for  thermofluid  structure  inter¬ 
action  applications  via  hydrodynamic  and/or  thermal  coupling 
mechanisms.  This  latest  development  has,  for  example,  facilitated  the  use 
of  coupled  discrete  element  method  (DEM)/LBM  models  that  can  be  used 
in  a  large  variety  of  applications  of  interest,  including  those  relating  to:  po¬ 
rous  media  flows  (Han  and  Cundall  1997)  (including  fines  migration  appli¬ 
cations),  geologic  deformation  (shear)  band  permeability  evolution  (Sun  et 
al.  2013),  gas-fluidized  beds  (Third  and  Miieller  2013),  and  liquid  phase 
sintering  (Varnick  et  al.  2013). 

Despite  the  increased  progress  and  stated  advantages  of  the  LBM  over 
conventional  approaches,  the  total  number  of  publications  resulting  from 
LBM  research  is  substantially  in  the  minority,  compared  to  more  conven¬ 
tional  approaches.  This  paucity  in  LBM  research  is  likely  due  to  its  relative 
novelty,  but  it  may  also  be  attributable  to  various  difficulties  relating  to 
implementation  issues,  including:  the  indirect  method  of  applying  bound¬ 
ary  conditions  (i.e.,  prescribing  distribution  functions  in  lieu  of  explicitly 
providing  Dirichlet  or  Neumann  conditions  involving  precise  macroscopic 
variables),  the  need  to  prescribe  aspect  ratio  and  Reynolds  number  (Re) 
equivalence  when  making  model  comparisons,  the  use  of  lattice  units  as 
compared  with  physical  (Re)  units,  or  the  numerical  instabilities  that  can 
arise  for  violations  of  compressibility  or  laminar  restraints.  Additional  lia¬ 
bilities  of  the  LBM  include  (Mohamad  2011): 

•  The  difficulty  in  simulating  high  Mach  number  (Ma)  flows  (in  excess  of 
the  incompressible  limit  -  e.g.,  Ma  >  0.3). 
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•  Local  violations  of  the  Kn.  which  are  particularly  likely  for  flows  with 
strong  velocity  gradients. 

•  The  limitations  associated  with  uniform  square  grids.  As  opposed  to 
conventional  FVM,  FEM,  or  EDM  approaches,  the  LBM  is  primarily  re¬ 
stricted  to  square  lattices.  Refinement  meshing  and/or  the  application 
of  curved  grids  have  proven  difficult  to  implement. 

•  The  LBM  is  subject  to  numerical  instabilities  particularly  when  the  lat¬ 
tice  viscosity  (v)  becomes  too  small. 

1.2  Objective 

While  there  is  currently  no  debate  concerning  the  theoretical  equivalence 
of  the  LBM  and  NS  for  low  Ma  flows  with  negligible  density  fluctuation 
(the  Chapman-Enskog  expansion  conclusively  proves  this,  as  discussed  in 
Chapter  4);  in  practice  however,  negligible  density  fluctuations  can  be  dif¬ 
ficult  to  achieve.  This  difficulty  is  particularly  true  for  imposed  pressure 
boundary  conditions.  In  fact,  regardless  of  the  imposed  boundary  condi¬ 
tions,  the  spatial  density  variation  is  never  completely  zero  in  LBM  simula¬ 
tions.  Other  potential  problems  underlying  the  traditional  LBM  are 
associated  with  local  stability  issues,  particularly  those  brought  on  by  local 
flow  gradients. 

In  light  of  the  aforementioned  concerns,  this  report  documents  a  compari¬ 
son/validation  effort  accompanying  the  development  of  a  standard  LBM 
solver  complemented  with  the  IMB  method.  Eor  validation  purposes,  con¬ 
ventional  NS  methods  are  utilized.  In  particular,  comparison  studies  will 
involve  four  standard  benchmark  cases  including,  the  flow  through  a  rec¬ 
tangular  channel,  the  flow  inside  a  lid  driven  cavity,  the  flow  over  a  back- 
step,  and  the  flow  over  a  stationary  circular  cylinder.  Each  case  will  be 
simulated  over  a  variety  of  different  input  conditions  and  associated 
Reynolds  numbers. 

1.3  Approach 

This  report  first  presents  the  incompressible,  laminar,  NS  equations 
(including  the  continuity  equation).  The  need  to  include  the  Poisson 
equation  is  demonstrated  for  purposes  of  closure.  Eor  illustration 
purposes,  a  simple  one-dimensional  convection-diffusion  example  is 
presented  and  discretized  in  accordance  with  the  EVM  of  solution.  Second, 
the  theoretical  background  for  the  LBM  is  presented  with  origins  derived 
from  the  standard  Boltzmann  equation.  In  lieu  of  the  complex 
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integrodifferential  collision  operator  of  the  standard  Boltzmann  equation, 
the  standard  BGK  approximation  is  used  (Higuera  and  Jimenz  1989; 
Koelman  1991).  Third,  assuming  negligible  density  variation,  the  analytical 
equivalence  between  the  LBM  and  NS  is  shown  using  the  Chapman- 
Enskog  expansion  method.  Finally,  the  aforementioned  benchmark 
comparison  studies  are  described  and  analyzed. 
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2  Navier-Stokes  Equations 

In  general,  the  state  of  most  fluids  can  be  described  by  using  the  following 
macroscopic  variables: 

u,  fluid  velocity 

p,  fluid  density 

p,  pressure 

T,  temperature 

E,  energy 

The  applicable  governing  equations  for  these  variables  can  be  derived  from 
statistical  mechanics  using  the  microscopic  equations  of  motion.  Alterna¬ 
tively,  they  can  be  obtained  through  conservation  laws.  In  divergence 
form,  the  conservation  of  mass  (Eq.  i),  momentum  (Eq.  2),  and  energy 
(Eq.  3)  can  be  expressed  as: 


Mass: 

|+v-(pn)  =  o 

(1) 

Momentum: 

—  ipu)  +  V  ■  {puu  +  pl  —  t)  =  pg 

(2) 

Energy: 

+  V  ■  (jipE  +  p)u  —  kVT  —  (u  ■  =  p(q  +  g  -  u) 

(3) 

where: 


V  =  gradient  operator 
g  =  acceleration  due  to  gravity 


q  =  internal  heat  source 


K  =  thermal  conductivity 
/  =  identity  matrix 


T  =  deviatoric  stress  in  the  fluid. 


Eor  a  Newtonian  fluid,  t  has  the  following  form  (Eq.  4): 


T  =  —p'piy  ■  u)I  +  2ppS 


(4) 
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Where  [i  and  pc'  represent  the  dynamic  (shear)  and  bulk  viscosity  respec¬ 
tively,  and  S'  (the  strain  rate  tensor)  is  defined  by  Eq.  5: 

5  =  1/2  (yu  +  (Vu)^)  (5) 

As  written,  for  an  n-dimensional  space,  Eq.  1-3  represent  a  system  of 
n  +  2  equations  for  n  +  4  unknowns.  To  close  the  system,  two  additional 
equations  of  state  are  required:  one  relating  the  pressure  and  temperature 
and  the  other  temperature  and  energy.  Eor  example,  for  a  perfect  gas, 
these  equations  are:  p  =  pRT  and  E  =  C^T,  respectively. 

2.1  Isothermal  and  incompressible  flows 

Eor  an  isothermal,  incompressible  flow,  the  effects  of  the  temperature  are 
neglected,  and  the  density  is  assumed  constant  (p  =  po).  The  conservation 
law  for  the  mass  is  simplified  and  states  that  the  velocity  field  is  diver¬ 
gence  free  (solenoidal),  shown  by  Eq.  6: 

V-u  =  0  (6) 

The  conservation  law  for  the  momentum  is  also  simplified  and  leads  to  the 
incompressible  NS  equations: 

dfU  +  (u  ■  V)u  =  —  ^Vp  +  vV^ii  (7) 

To  find  the  value  of  pressure,  p,  one  takes  the  divergence  of  Eq.  7.  Making 
use  of  Eq.  6,  this  yields  the  Poisson  equation  (Eq.  8): 

V2p  =  -po(Vu):(Vuf  (8) 

Eor  implementation  purposes,  Eq.  8  is  solved  via  an  iterative  procedure  at 
each  time  step.  The  correct  value  of  pressure  is  obtained  in  a  way  that  en¬ 
sures  the  velocity  field  remains  divergence-free. 

2.2  Finite  volume  discretization 

While  there  are  many  methods  that  can  be  used  to  discretize  the  above 
conservation  equations,  for  this  work  we  implement  the  EVM  (Versteeg 
and  Malalasekera  1995).  Eor  illustration  purposes,  the  EVM  discretization 
pertaining  to  a  steady-state  flow  involving  one-dimensional  convection 
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and  diffusion  of  a  quantity  {cp),  is  demonstrated  for  the  control  volume 
(Figure  i). 


Figure  1.  One-dimensional  control  volume  used  in  the  FVM. 
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For  a  diffusion  coefficient  (F),  Eq.  6  and  Eq.  7  become: 


(9) 


—  (ucp)  = 

dx  dx  \  dxj 


(10) 


Integration  of  Eq.  9  and  Eq.  10  over  the  control  volume  gives  us  Eq.  11  and 
Eq.  12: 


(uA(p)e  -  {uA(p)^  =  0 


(11) 


-  (UAV)„  =  (pa  -  (pa  (12) 


where: 

e  and  w  represent  directional  indices  as  shown  in  Eigure  1,  and  A  is  the 
control  volume  cross-sectional  area. 

As  customary,  we  define  the  variables  (Fg,  and  Dg,  to  correspond  to 
convection  and  diffusion  in  the  east  and  west  directions  in  Eq.  13  and 
Eq.  14,  respectively: 


(^)e> 


fw  (^)w 


(13) 
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SXy\rp 


De 


SXpE 


(14) 


Here  the  indices  W,  P,  and  E  correspond  to  the  respective  node  locations 
shown  in  Figure  i. 


Substituting  Eq.  13  and  Eq.  14  into  Eq.  11  and  Eq.  12,  and  assuming  that  all 
cross-sectional  areas  are  equivalent  (i.e.,  =  Ag  =  A),  we  have  Eq.  15  and 

Eq  16: 


Fe-Fy,  =  0  (15) 

Fe(pe  -  FwCpw  =  DgicpE  -  (pp)  -  Dyy((pp  -  cpw)  (16) 

Einally,  the  convection  terms  on  the  left-hand  side  of  Eq.  i6  may  be  treat¬ 
ed  with  an  upwind  differencing  scheme,  which  takes  into  account  the  di¬ 
rection  of  the  flow.  Eor  positive  flows  (streaming  from  left  to  right),  we 
have  Eq.  17: 

<Pw  —  (Pw  q)g  =  (pp  for-.  Uyy  >  0,  and  Ug  >  0  (17) 

and  for  negative  flows  (right  to  left),  we  have  Eq.  18: 

(Pyy  =  q)p,  and  cpg  =  cpp  for-.  Uyy  <  0,and  Ug  <  0  (18) 

This  completes  the  discretization  process  for  this  simple,  one-dimensional, 
convection-diffusion  case. 
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3  Lattice  Boltzmann  Equation 

As  its  name  suggests,  the  Lattice  Boltzmann  (LB)  equation  may  be  regard¬ 
ed  as  a  limited  version  of  the  generalized  Boltzmann  equation.  The  Boltz¬ 
mann  equation  describes  the  phenomenological  properties  of  particulate 
systems  via  the  use  of  the  phase  space  distribution  function,  f(x,  c,  t).  The 
phase  space  distribution  function  represents  the  number  of  particles  at 
time  t,  positioned  between  x  and  x  +  dx,  with  velocities  between  c  and 
c  +  dc.  As  shown  in  Figure  2,  if  we  impose  an  external  force  (F)  to  a  parti¬ 
cle  of  unit  mass,  this  force  will  result  in  a  change  in  position  and  velocity, 
equivalent  to:  x  +  cdt  and  c  +  Fdt,  respectively. 

Figure  2.  Position  and  veiocity  vectors  resuiting  from  an  externai  force. 


y 


For  the  idealized  case  in  which  there  are  no  collisions  between  particles, 
we  may  equate  the  particle  distribution  at  time  (f)  with  that  of  time 
(t  +  dt),  to  get  Eq.  19: 

f{x  +  cdt,  c  +  Fdt,  t  +  dt)dxdc  —  f{x,  c,  t)  =  0  (19) 

If  collisions  are  included,  it  is  clear  that  some  particles,  initially  at  (x,  c,  t), 
will  not  arrive  at  (x  +  cdt,  c  +  Fdt,  t  +  dt)  simply  because  they  have  been 
diverted  from  their  original  path.  Analytically,  this  may  be  represented  as 
Eq.  20: 

/(x  +  cdt,  c  +  Fdt,  t  +  dt)dxdc  —  /(x,  c,  t)drdc  =  n(f)dxdcdt  (20) 
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where  we  have  used  n(/)  to  represent  the  collision  operator  (physically 
representing  the  difference  from  the  idealized  streaming  case).  Dividing 
Eq.  20  by  dxdcdt  and  taking  the  limit  as  dt  ^  0  yields  Eq.  2i: 

f  =  m  (21) 

Since/is  a  function  of  x,  c,  and  t,  the  total  rate  of  change  of/may  be  writ¬ 
ten  as  Eq.  22: 


df  =  (yxHdx  +  (Vc/)dc  +  ^dt  (22) 

Dividing  by  dt  and  applying  Eq.  21,  we  obtain  the  familiar  Boltzmann 
equation  shown  as  Eq.  23  or  Eq.  24: 

^  ^  +  (V,/)c  +  /m  =  n(0  (23) 


OR: 

^  +  (?'V*  +  ^-y  +  a,)/  =  n(f)  (24) 

Where  Vx  (for  the  two  dimensional,  Cartesian  coordinate  system  shown  in 
Figure 2) is (£,^)  andV,is(|;,^). 

The  LB  equation  is  easily  derived  from  Eq.  24.  By  neglecting  external  forc¬ 
es  (F),  normalizing  the  mass  to  unity,  and  limiting  the  velocity  (c)  to  a  dis¬ 
crete  set  of  vectors  (in  conformance  with  a  prescribed  lattice),  the 
distribution  function  can  be  written  as  /j  (t,  x)  with  no  loss  of  information 
(Mohamad  2011).  With  time  discretization,  we  have:  fldt  fli(r,  t),  and 
the  LB  equation  can  be  written  as  Eq.  25: 

fi{t  +  5t,x  +  CiSt)  -  fi{x,t)  =  nj(x,t)  (25) 

This  equation  has  been  written  in  lattice  units,  such  that  the  time  interval 
from  one  iteration  to  the  next  is  and  the  spacing  between  two  adjacent 
lattice  nodes  is  Sx,  where  Sx  —  CiSf-. 

As  stated  previously,  one  of  the  main  problems  in  solving  the  Boltzmann 
equation  (Eq.  24)  is  the  complicated  nature  of  the  integrodifferential  colli- 
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sion  term  (Allen  2006)  The  current  practice  within  the  context  of  the  LBM 
is  to  use  the  BGK  relaxation  approximation  for  the  collision  term  (Higuera 
and  Jimenez  1989;  Koelman  1991).  The  BGK  approximation  involves  re¬ 
laxation  dynamics  (with  relaxation  parameter  co)  towards  local  equilibri¬ 
um,  to  produce  Eq.  26: 

=  (26) 

As  shown  in  Figure  3,  the  D2Q9  particle  velocity  model  is  used  in  this 
work  (Qian  et  al.  1992).  As  shown,  the  particle  velocities  consist  of  a  rest 
particle  (cq)  and  eight  moving  particles  (c^,  C2, ...  Cg). 

Figure  3.  Lattice  veiocities  corresponding  to  the  D2Q9  configuration. 


The  local  equilibrium  distribution  function  (/j^^)is  obtained  from  a  series 
approximation  (see  Appendix)  and  defined  (up  to  second  order)  as  Eq.  27: 


/7eq  _  *. ..  ii  I  ^ia'^a  , 


The  speed  of  sound  constant  (c^)  and  the  lattice  weighting  coefficients  (wj) 
are  determined  such  that  a  system  of  moment  summations  are  satisfied. 
For  a  D2Q9  lattice  configuration  (Qian  et  al.  1992),  these  include  Eq.  28- 

31): 


(28) 


Zifr-p 

=  P'^a 


(29) 
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lif^CiaCip  =  \pCs'^Sap  +  pU^Up  (30) 

Si/t  ^ia^ip^iy  ~  3^^^  {^ap^y  ”1”  ^ps^a  ”1”  ^yaP^p')  (31) 

Where  the  subscript  indices  (a,  p,...)  represent  the  spatial  components  of 
the  vectors,  Cg  is  the  lattice  speed  of  sound,  and  5  is  the  Kronecker  delta. 

Satisfying  Eq.  28-31,  a  value  of  c|  =  1/3  is  often  used,  with  weighting  coef¬ 
ficients  as  shown  in  Eq.  32: 

Wj  =  ^(i  =  0),  Wj  =  i(i  =  1,2, 3, 4),  and  Wj  =  ;^(i  =  5,6,7,8)  (32) 

9  9  36 

The  macroscopic  density  (p)  and  velocity  (Uq.)  are  calculated  from  moment 
summations  using  the  distribution  function,  as  shown  in  Eq.  33  and  Eq. 

34,  respectively: 


P  =  lift 


(33) 


P^a  ~  Yii  fi^ia 


(34) 
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4  Lattice  Boltzmann/Navier-Stokes 

Equivalence  using  the  Chapman-Enskog 
Expansion 

Under  the  assumption  of  negligible  density  fluctuations  (e.g.,  incompress¬ 
ible  flow),  the  Chapman-Enskog  expansion  (Chapman  and  Cowling  1970) 
can  be  used  to  show  LBM  equivalence  to  the  NS  equations.  To  summarize, 
the  basic  idea  behind  the  Chapman-Enskog  expansion  is  to  separate  the 
distribution  function  into  multiple  scales  with  respect  to  the  order  of 
Knudsen  number  (e).  Physical  properties  of  the  macroscopic  variables, 
such  as  the  density  and  momentum,  are  automatically  separated  out  from 
the  different  scales. 

In  detail,  the  process  begins  first  with  the  distribution  function  written  as 
an  asymptotic  series  expanded  over  the  Knudsen  number  (e),  as  in  Eq.  35: 

/i  =/® +  +  +  -  (35) 

The  functions  (//^^)  are  defined  in  such  a  way  as  to  progressively  and  in¬ 
dependently  tend  toward  zero  with  increasing  powers  of  epsilon  (which 
necessarily  assumes  that  epsilon  is  a  very  small  quantity  consistent  with 
the  continuum  assumption).  Since  the  main  objective  of  the  Chapman- 
Enskog  expansion  is  to  provide  a  consistent  definition  for  ,  preliminar¬ 
ily  this  means  satisfying  Eq.  33  and  Eq.  34  as  Eq.  36  and  Eq.  37  by  requir¬ 
ing  that  the  first  two  moments  of  the  zeroth  approximation  reproduce 
macroscopic  density  and  velocity.  The  corresponding  moments  of  the 
higher-order  terms  are  set  to  zero: 


(36) 

Si/,®  =0; 

for  k>0 

(37) 

Using  the  BGK  approximation,  the  LB  equation  may  be  written  as  Eq.  38: 
/j(t  +  At,x  +  CiSt)  -fiit.x)  =  -i[/j(t,x)  (i  =  0  ...N) 


(38) 
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Applying  a  Taylor  series  expansion  (t,  x)  to  the  first  term  of  Eq.  38  yields 
Eq.  39: 


f,(t  +  S,J  +  cA)  =/<  +  «.  (s  +  A)  /<  +  +  om  (39) 


Where  fi(t,  x)  is  simplified  to  /j  for  convenience,  and: 

fd  a  /  32  32  f 

\dt  ^  dxj  J ‘  “  V3t2  ^  at3x„  ^  j  ^ ‘ 


(40) 


Substituting  Eq.  39  into  Eq.  38  and  dividing  by  Sf-  yields  Eq.  41: 

{h  +  +  7  (^  +  ^ 

Since  e  and  5^-  are  approximately  the  same  magnitude  (for  Cg  =  0(1)),  the 
distribution  function  /j  can  be  expanded  in  the  asymptotic  series  shown  in 
Eq.  42: 


ft  =  ft"  +  s,ft‘‘  +  +  - 


O)  X  2  /r(2) 


(42) 


Inserting  Eq.  4)  into  Eq.  41,  we  have  Eq.  43: 


dXr 


St(d 
^  2  (dt  ^ 


dx, 


)(/i'‘'  +  «t/™  +  «t"/i® +  ■•■) 

)  (/)'■' +«t/i“  +  '5tV“  +  -) 


+  5)7(^'’  +  “  •'■°’)  = 

Organizing  terms  according  to  the  order  of  Sf-  we  have  Eq.  44: 

d  \  1 


d 


^eq 


+  7fi 


+5f 


V3t  ^2V3t 


(2) 


=  0(5,^) 


(43) 


(44) 
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To  obtain  the  Euler  expression  (the  inviscid  form  of  NS),  the  first  term  of 
Eq.  44  is  assumed  zero,  creating  Eq.  45: 


“  (45) 

Taking  the  summation  Ej  on  Eq.  (45)  and  using  Eq.  28,  Eq.  29,  and  Eq.  36, 
we  obtain  the  equation  for  the  conservation  of  mass  in  Eq.  46: 


=  HP +  5^  »  (46) 

Taking  the  summation  'Zi  on  Eq.  45  and  using  Eq.  29,  Eq.  30,  and 
Eq.  37,  we  obtain  Euler’s  equation  for  the  conservation  of  momentum  (Eq. 

48): 


=  J^ipy-a)  +  +  pUaUp)  =  0 


(47) 


Next  the  contribution  from  the  second  term  of  Eq.  (44)  is  evaluated.  Ap¬ 
plying  the  differential  operation  (-^  +  CjQ.  on  Eq.  (45)  we  have  Eq.  48: 

\CfL  OOCrj// 


(48) 
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Inserting  Eq.  48  into  the  second  term  of  Eq.  44  to  eliminate  the  squared 
differentiation  results  in  Eq.  49: 


Taking  the  summation  Ej  on  Eq.  49  and  using  Eqs.  36  and  37,  we  have  Eq. 
50: 


= [(1  -  h)  (tM'" + "'“) + “  <50) 

Taking  the  summation  'Zi  CtaOn  Eq.  (49)  and  inserting  Eq.  37  yields  Eq.  51: 


=  At 


(1  ■  h)  (sZ  +  a^Z  S  ^ 


r(2) 

!■ 


(51) 


Erom  Eq.  45,  the  distribution  is  shown  in  Eq.  52: 


eq 


(52) 


Therefore  the  term  Zi  CjQ.Cj^  in  Eq.  51  is  transformed  to  Eq.  53,  using 
Eq.  30  and  Eq.  31: 


+  c 


dXyl'^ 


f(eq) 

Ti 


=  — T 


d_(l 


—  — T 


Y^P  +  +  ^Py'^a  + 


(53) 
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Applying  the  product  rule  to  the  term  ^  {^pu^up)  of  Eq.  53  yields  Eq.  54: 

^  {pUaUp)  =  Up  ipUa)  +  ^  {pUp)  +  U^Up  (54) 

Using  Eq.  46  and  Eq.  47  for  and  ^  ipu^),  respectively,  we  have  Eq.  55: 


— T 


) 


(55) 


Summing  Eq.  47  and  Eq.  51  with  Eq.  55  and  neglecting  ^  (^pUaUpUy)  as  a 

small  quantity  results  in  the  momentum  equation  (Eq.  56)  in  the  NS  equa¬ 
tions: 


d  d 


d 

dxp 


fdUg  ^  dup\ 
\dxp  dXa) 


=  0 


(56) 


Where  the  pressure  (P)  and  viscosity  (p)  are  expressed  as  Eq.  57  and  Eq. 
58,  respectively: 


=  \pc^ 

(57) 

(58) 
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5  Implementation  Issues 

As  shown  in  the  previous  section,  the  incompressible  NS  equations  can  be 
derived  from  the  LB  equation  through  the  Chapman-Enskog  procedure  if 
the  density  fluetuation  is  assumed  negligible.  In  practice  however,  this 
last  assumption  is  never  completely  true.  Indeed,  the  spatial  density  varia¬ 
tion  is  never  completely  zero  in  LBM  simulations.  This  is  particularly  true 
for  flow  applications  wherein  a  pressure  gradient  is  imposed  at  the  bound¬ 
aries.  Although  there  have  been  recent  efforts  to  reduce  or  completely 
eliminate  the  compressible  effect  in  the  LBM,  the  results  have  not  been  en¬ 
tirely  satisfactory,  particularly  for  unsteady  flows  (Frisch  et  al.  1987;  Alex¬ 
ander  et  al.  1992;  Zou  et  al.  1995). 

Despite  its  intrinsic  compressibility,  the  LBM  is  applicable  only  to  the  low 
Mach  number  (Ma  =  u/Cg)  flow  applications.  As  shown  in  the  previous 
section,  this  is  because  a  small  velocity  expansion  is  (implieitly)  used  with¬ 
in  the  Chapman-Enskog  expansion.  This  low  Mach  number  approximation 
is  equivalent  to  the  incompressible  limit  (i.e.,  Ma  <  0.3). 

To  correctly  simulate  incompressible  flow  in  practice,  one  must  ensure 
that  the  Mach  number,  and  the  density  variation,  Sp,  are  sufficiently  small, 
of  order  0(6)  and  0(6^),  respectively,  where  the  Knudsen  number  may  be 
approximated  as  e  =  (ct)/L.  Quantitatively,  this  can  be  assured  from  the 
expression  of  the  Mach  number  as  a  function  of  the  Reynolds  number. 
Since  the  fluid  viscosity  is  related  to  the  relaxation  frequency  (co  =  1/t), 
we  have  Eq.  59: 


v  =  1^(0) -0.5)  (59) 

The  Mach  number,  Ma,  can  be  obtained  by  dividing  both  sides  of  equation 
X  by  UL.  This  action  yields  Eq.  60: 

Ma  =  ^(o)  -  O.S)Re  (60) 

In  Eq.  60  L/Sy.  is  the  number  of  lattice  sites,  N,  in  the  direction  of  the 
characteristic  length  (normally  perpendicular  to  the  flow  direction),  and 
Re  is  the  Reynolds  number  (Re  =  UN/  v).  Since  Sy  is  normally  assigned  a 
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value  of  unity  (for  purposes  of  convenience)  to  keep  Ma  small  (Ma  <  0.3), 
then  0),  N,  and  U  should  be  selected  accordingly.  In  practice,  U  is  generally 
prescribed  values  on  the  order  of  0.1.  Since  Reynolds  number  similarity 
must  also  be  enforced  between  the  macroscale  and  lattice  environments, 
this  enforcement  places  an  additional  constraint  on  the  selection  of  o)  and 
N. 

In  addition  to  the  aforementioned  constraints,  additional  stability  issues 
may  arise  due  to  the  presence  of  parasitic  oscillations  in  the  vicinity  of 
sharp  gradients  in  the  flow.  These  sharp  gradients  may  be  caused  by  shock 
waves  travelling  through  the  fluid  or  by  excessively  thin  shear  layers. 

While  conventional  methods  allow  for  a  certain  amount  of  numerical,  arti¬ 
ficial  diffusion  to  mitigate  impending  instabilities  of  this  sort,  this  mitiga¬ 
tion  is  not  the  case  for  the  traditional  LBM.  Indeed,  as  shown  in  Eq.  26, 
the  traditional  BGK  collision  approximation  involves  only  a  single,  con¬ 
stant-valued  relaxation  time.  While  several  novel  “dissipative”  methods 
have  recently  attempted  to  address  this  stabilization  issue,  including  colli¬ 
sion  operations  involving  multiple  relaxation  times  (MRTs)  (d’Humieres 
1994;  d’Humieres  et  al.  2002),  these  methods  have  demonstrated  only 
limited  success.  This  lack  of  success  is  largely  due  to  the  fact  that  although 
the  dissipative  methods  can  provide  substantial  (or  complete)  oscillatory 
suppression,  accuracy  is  compromised.  However,  due  to  the  low  Mach 
number  used  in  this  work  in  incompressible  flow  restraints,  this  type  of 
stability  issue  is  of  negligible  concern. 
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6  Lattice  Boltzmann  and  Navier-Stokes 
Benchmark  Applications 


The  following  benchmark  applications  were  selected  for  their  ubiquity  in 
the  scientific  and  engineering  literature,  relative  ease  to  set-up,  and  com¬ 
putational  efficiency.  All  of  the  simulations  were  run  in  two  dimensions 
and  in  serial  using  a  Dell  Precision  490  desktop  machine  with  two  Dual- 
Core  Xeon  64  bit  processors.  Post  processing  was  facilitated  using  VTK*- 
formatted  output  files  and  Paraviewt  (Henderson  2007). 


The  Navier-Stokes  simulations  were  conducted  using  the  finite  volume 
solver,  OpenFOAM  version  2.1.  *  In  particular,  the  incompressible  laminar 
flow  equations  (Eq.  6-8)  were  solved  in  conjunction  with  the  pressure  im¬ 
plicit  with  splitting  operator  (PISO)  algorithm  (Issa  1986).  Discretization 
of  diffusion  and  convective  terms  was  carried  out  using  standard  second- 
order  central  differencing  and  upwinding  schemes,  respectively.  Since  the 
OpenFOAM  code  is  inherently  transient,  steady-state  conditions  were  ob¬ 
tained  over  extended  time  periods,  as  quantified  by  residual  values  of  less 
than  10-6.  Specifically,  the  steady-state  condition  was  checked  by  calculat¬ 
ing  the  residual  (Res)  of  the  x-component  of  velocity  at  each  time  step,  in 
accordance  with  the  following  equation: 


Res  = 


vrr  yLX  f 
Zjy=i  Zjx=i  y 


\u(x,y)t+i-u(x,y) 


\u(x,y)t 


(61) 


The  LBM  code  routines  were  developed  in-house  by  ERDC  and  pro¬ 
grammed  in  FORTRAN  90  using  much  of  the  prescribed  content  in  Chap¬ 
ter  3.  Collaborations  with  Mississippi  State  University’s  Center  for 


*  The  Visualization  Toolkit  (VTK)  is  an  open-source,  freely  available  software  system  for  3D  computer 
graphics,  image  processing  and  visualization. 

t  Paraview  is  an  open-source  multiple-platform  application  for  interactive,  scientific  visualization,  which 
began  in  2000  as  a  collaborative  effort  between  Kitware  and  Los  Alamos  National  Laboratory. 

OpenFOAM  is  a  free,  open-source  computational  fluid  dynamics  software  developed  by  OpenCFD  Ltd. 
of  the  ESI  Group  and  distributed  by  the  Open  Foam  Foundation.  More  information  is  available  at 
openfoam.org/version  2.1.1. 
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Advanced  Vehicular  Systems  (CAVS)  were  of  significant  benefit  in  the  code 
development  process.  A  two-dimensional  square  lattice  with  nine  lattice 
velocities  (D2Q9)  was  used  (Qian  et  al.  1992).  Boundary  conditions  (in 
particular  the  bounce  back,  velocity  inlet/outlet  conditions,  and  pressure 
conditions),  were  implemented  in  accordance  with  those  prescribed  by 
Zou  and  He  (1997).  In  addition  to  the  standard  LBM,  the  implementation 
of  the  IMB  method  (Owen  et  al.  2011)  was  also  incorporated  for  the  case  of 
the  flow  over  a  circular  cylinder.  Implemented  in  accordance  with  Owen  et 
al.  (2011),  this  method  better  represents  the  intersection  of  the  fluid  with 
the  curved  cylinder  walls.  While  the  exact  details  are  omitted  in  this  report 
for  reasons  of  scope,  Eq.  62-65  outline  the  basic  procedure. 

The  standard  LBM  equation  (Eq.  38)  is  modified  according  to  Eq.  62: 

fi(x  +  BiSt,  t  + St)-  fi(x,  t)  =  (1  -  Bn)  t)  -  fi(x,  t))  j  +  BnD.f  (62) 

Where  Bn  is  an  empirically  derived  weighting  function  of  the  area  (or  vol¬ 
ume)  fraction  (f)  of  lattice  cell  (i)  that  is  occupied  by  the  cylinder,  the  re¬ 
sult  is  Eq.  63: 


Bnis)  = 


£(t-1/2) 
(1-£)  +  (t-1/2) 


(63) 


Where  Hf  is  a  collision  operator  representing  the  change  of  momentum  due 
to  collisions  with  the  cylinder,  we  have  Eq.  64: 

t)  -  v^)]  -  [/i(x,  t)  -  Vp)]  (64) 

Where  i'  represents  the  direction  opposite  Ij,  and  Vp  is  the  velocity  at  the 
surface  of  the  cylinder  at  point  p  (located  at  close  proximity  to  lattice  node 
i),  we  have  Eq.  65: 


^cm  L  aj  X  Tpjf-ni  (65) 

Where  Vcm  is  the  velocity  of  the  cylinder’s  center  of  mass,  cj  is  the  angular 
velocity,  and  rp/cm  is  the  position  vector  atp  relative  to  the  center  of  mass 

of  the  cylinder. 
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6.1  Rectangular  channel  flow 

6.1.1  Problem  description  and  setup 

The  incompressible  flow  through  a  rectangular  channel  represents  a  com¬ 
mon  computational  fluid  dynamics  (CFD)  benchmark  problem.  The  flow  is 
typically  driven  by  either  a  constant  inlet  velocity  condition  or  by  means  of 
a  pressure  gradient.  In  this  particular  application,  water  flows  between 
two  parallel  plates  of  length  (L  =  0.5  m),  height  (H  —  0.02  m)  and  of  aspect 
ratio,  AR  =  25.  Prescribed  velocity  /pressure  gradient  conditions  are  ap¬ 
plied  at  the  inlet  and  outlet  boundaries,  while  no-slip  (zero  velocity)  condi¬ 
tions  are  applied  along  the  walls. 

For  the  velocity  inlet  condition,  three  different  Reynolds  numbers  are  in¬ 
vestigated  using  both  the  LB  and  NS  methods.  These  cases  correspond  to 
Re=200,  Re=400,  and  Re=8oo.  The  pressure  gradient  condition  is  also 
examined  using  both  methods,  but  only  for  Re  =  10.  The  corresponding 
input  parameters  are  shown  in  Table  1,  with  a  constant  lattice/grid  resolu¬ 
tion  of  1000x40  maintained  throughout. 


Table  1.  Simulation  parameters  used  for  the  rectangular  channel  flow  problem. 


Case 

Inlet  Vel. 

Nav-Stokes 

[m/s] 

Inlet  Vel.  (f/^) 
LB  [l/t] 

Kin.  Vise,  (v), 

Nav-Stokes 

[m2/s] 

Kin.  Vise. 
(V),  LB 
[P/t] 

AP 

Re  =  10 

— 

— 

l.OE-6 

0.16667 

l.OE-5 

Re  =  200 

0.01 

0.3 

l.OE-6 

0.01 

— 

Re  =  400 

0.02 

0.2 

l.OE-6 

0.02 

— 

Re  =  800 

0.04 

0.1 

l.OE-6 

0.01 

— 

6.1.2 

Results 

The  velocity  profiles  comparing  the  results  of  the  LBM  and  NS  method 
when  using  various  velocity  inlet  conditions  are  shown  in  Figure  4.  As  ex¬ 
pected,  the  velocity  distribution  along  the  downstream  direction  progres¬ 
sively  evolves  into  a  characteristic  parabolic  profile  as  the  flow  reaches  a 
state  of  steady  equilibrium.  Also  as  expected,  increasing  the  Reynolds 
number  results  in  higher  velocity  profiles  along  each  of  the  downstream 
positions  (X/H  =  o.iL,  X/H=o.2L,  andX/H=o.5L).  Comparisons  between 
the  LBM  and  NS  methods  show  excellent  agreement. 
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The  velocity  profiles  comparing  the  results  of  the  LBM  and  NS  method 
when  using  a  pressure  gradient  boundary  condition  are  shown  in  Figure  5. 
Unlike  Figure  4,  the  imposed  pressure  gradient  results  in  a  uniform  veloci¬ 
ty  distribution  along  each  of  the  downstream  sampling  locations.  Like  Fig¬ 
ure  4,  the  comparisons  between  the  LBM  and  NS  methods  again  show 
excellent  agreement. 


Finally,  Figure  6  shows  a  representative  time  history  of  the  convergence 
error  in  accordance  with  Eq.  61.  As  indicated,  the  residual  for  x- 
component  of  velocity  becomes  fully  convergent  (reaches  a  steady-state 
condition)  at  approximately  2,000  iterations. 

Figure  4.  Velocity  profiles  comparing  the  results  of  the  LBM  and  NS  corresponding  to 
the  flowthrough  a  rectangular  channel  (a)  ffe=  200,  (b)  ffe=  400,  and  (c)  /?e=  800). 
A  velocity  inlet  condition  is  used.  As  indicated,  the  velocities  were  sampled  at  three 
downstream  locations  corresponding  to  X/L  =  0.1,  X/L  =  0.2  and  X/L  =  0.5. 
(Lattice/grid  resolution  =  1000x40). 


(a)  Re  =  200 
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(b)  Re  =  400 


(c)  Re  =  800 
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Figure  5.  Velocity  profiles  comparing  the  results  of  the  LBM  and  NS  corresponding  to 
the  flow  through  a  rectangular  channel  {Re  =  200,  Re  =  400,  and  Re  =  800).  A 
pressure  gradient  condition  is  used.  As  indicated,  the  velocities  were  sampled  at 
three  downstream  locations  corresponding  to  VL  =  0.1,  VL  =  0.2  and  VL  =  0.5. 
(Lattice/grid  resolution  =  1,000x40). 


Figure  6.  Representative  velocity  convergence  error  per  Eq.  61  for  the  flow  through  a 
rectangular  channel.  As  shown,  approximately  2,000  LBM  iterations  were  required  for 

fully  steady-state  conditions. 
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6.2  Flow  through  a  lid-driven  cavity 

6.2.1  Problem  description  and  setup 

The  lid-driven  cavity  represents  another  common  CFD  benchmark  prob¬ 
lem.  The  moving  lid  creates  a  strong  vortex  in  the  center  of  the  domain, 
and  a  series  of  weaker,  secondary  vortices  in  the  lower  left  and  right  cor¬ 
ners  are  typical  of  successively  higher  Reynolds  numbers.  In  this  particular 
application,  a  square  cavity  with  unit  aspect  ratio  {L  =  H  -  2m)  is  filled 
with  engine  oil  at  15°C  (v  =  1.2E-3  m^/s).  The  lid  is  set  in  motion  (provid¬ 
ing  a  fluid  momentum  source),  and  zero  velocity  (no-slip)  boundary  condi¬ 
tions  are  maintained  along  the  remaining  walls. 

For  this  problem,  three  different  Reynolds  numbers  are  investigated 
(Re=ioo,  Re=i,ooo,  andRe=3,ooo),  and  the  corresponding  input  param¬ 
eters  are  shown  in  Table  2  for  both  the  LBM  and  NS  methods.  A  constant 
lattice/grid  resolution  of  100x100  is  maintained  for  each  case. 


Table  2.  Simulation  parameters  used  for  lid-driven  cavity  problem. 


Case 

LidVel.  (Uud). 
Nav-Stokes  [m/s] 

Inlet  Vel. 

{Uiat)  LB  [l/i] 

Kin.  Vise,  (v),  Nav- 
Stokes  [mys] 

Kin.  Vise,  (v),  LB 
[lyt] 

Re  =  100 

0.6 

0.1 

1.2E-3 

0.1 

Re  =  1000 

6 

0.1 

1.2E-3 

0.01 

Re  =  3000 

18 

0.01 

1.2E-3 

0.01 

6.2.2  Results 

Comparisons  of  the  velocity  contours  are  shown  in  Figure  7.  As  indicated 
for  the  prescribed  simulated  conditions,  both  methods  show  remarkable 
qualitative  agreement  and  reveal  the  presence  of  a  large  central  vortex 
which  tends  to  migrate  to  the  upper  right  corner  for  decreasing  Reynolds 
numbers.  Also  observed  is  the  presence  of  two  secondary  (counter¬ 
rotating)  vortices  located  near  the  lower  left  and  right  of  the  domain.  As 
shown,  these  tend  to  increase  in  size  with  increasing  Reynolds  number. 
For  Re  =  3,000,  an  additional  secondary  vortex  develops  along  the  upper 
left  side. 

Quantitative  comparisons  of  velocity  results  between  the  LBM  and  N-S 
method  are  shown  in  Figure  8.  The  results  correspond  to  the  three  afore¬ 
mentioned  Reynolds  numbers,  and  mesh/lattice  resolutions  of  100x100. 
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As  expected,  in  the  vertical  direction  (Figure  8(a)),  the  general  x-velocity 
profiles  show  maximum  positive  velocities  near  Y/H  =i,  a  direction  change 
at  approximately  Y/H=o.5,  and  maximum  negative  x-velocities  for 
0.1  <Y/H<  0.2.  Correspondingly,  in  the  horizontal  direction  (Figure 
8(b))  the  y  component  of  velocity  shows  maximum  positive  velocities  for 
0.1  <  X/L  <  0.2,  a  change  in  direction  near  A/L  =  0.5,  and  maximum  neg¬ 
ative  velocities  at  approximately  0.89  <  XfL  <  0.95.  Comparisons  between 
the  two  methods,  along  the  vertical  and  horizontal  centerlines,  are  in  gen¬ 
eral  excellent  and  reveal  minimal  differences,  particularly  for  i?e  =  100, 
and  i?e=i,ooo.  For  Re=3,ooo  however,  some  disparity  is  observed  particu¬ 
larly  along  the  upper  and  lower  limits  of  Y/H.  This  is  likely  the  result  of  the 
increased  compressibility  effects  associated  with  larger  Reynolds  numbers 
and  the  limitations  of  the  LBM  in  handling  these  cases  (as  described  in 
Chapter  5). 

Figure  7.  Velocity  contours  comparing  the  results  of  the  LBM  and  Navier-Stokes  for 
flow  through  a  lid-driven  cavity  {Re  =  100,  Re  =  1,000  and  Re  =  3,000).  As  shown, 
the  presence  of  a  large  central  vortex  is  observed  for  each  case, 
as  well  as  several  secondary  vortices  (lattice/grid  resolutionilOOxlOO). 


(a)  NS;Re=ioo 


(b)  LBM;  Re=ioo 
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I  I  I  I  I  I  I  I 


(c)  NS;  i?e=l  ,000  (d)  LBM;  7?e=l  ,000 


(e)  NS;i?e=3,000 


(f)  LBM:  i?e=3,000 


Y/H 
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Figure  8.  Velocity  profiles  comparing  the  results  of  the  LBM  and  NS  for  flow  through 
a  lid-driven  cavity  {Re=  100,  Re=  1,000  and  /?e=  3,000).  As  indicated,  the  x- 
component  and  y-component  of  velocity  were  sampled  along  the  vertical  and 
horizontal  domain  centerlines,  respectively. 


(a)  Vertical  centerline 


(b)  Horizontal  centerline 


6.3  Rectangular  channel  flow  with  back-step 

6.3.1  Problem  description  and  setup 

A  third  CFD  benchmark  problem  investigated  in  this  report  involves  the 
laminar  incompressible  flow  over  a  back-step.  In  contrast  to  turbulent 
flows  (wherein  a  single  recirculation  zone  is  created  just  aft  of  the  step), 
laminar  flows  exhibit  various  recirculation  zones  occurring  downstream  of 
the  step.  Flow  separation  occurs  when  adverse  pressure  gradients  act  on 
the  fluid.  For  a  low-to-moderate  initial  Reynolds  number,  the  flrst  region 
of  separation  occurs  just  aft  of  the  step  along  the  bottom  wall.  With  in¬ 
creasing  Reynolds  number,  a  second  region  of  separation  occurs  along  the 
top  wall.  Further  increases  in  the  Reynolds  number  create  yet  a  third  sepa¬ 
ration  region  downstream  of  the  first  recirculation  region  along  the  bot¬ 
tom  wall.  In  fact,  recirculation  zones  will  continue  to  develop  downstream 
in  this  manner  as  long  as  the  Reynolds  number  continues  to  increase  and 
the  flow  remains  laminar.  Of  course  this  is  rarely  observed,  as  the  flow 
eventually  becomes  turbulent  (Jongebloed  2008). 

The  benchmark  problem  used  in  this  study  consists  of  a  channel  with  a 
backward-facing  step  at  the  entry  (Figure  9).  As  shown,  the  height  of  the 
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back-step  is  half  that  of  the  channel,  and  extends  one  channel  height  (H  = 
0.02  m)  before  it  opens  to  the  flow.  The  length  of  the  channel  {L=25H)  is 
long  enough  for  the  flow  to  become  fully  developed  (as  determined  from 
the  rectangular  channel  flow  discussion  in  Section  6.1).  The  flow  is  driven 
by  a  constant  inlet  velocity  condition  with  no-slip  (zero  velocity)  condi¬ 
tions  at  the  upper  and  lower  walls. 

For  this  problem,  comparisons  between  the  LB  and  NS  methods  are  con¬ 
ducted  for  three  different  Reynolds  numbers  (Re  =  200,  Re  =  400,  and  Re 
=  800).  An  additional  simulation  is  conducted  for  Re  =  1,000  (using  solely 
the  LBM)  to  illustrate  the  vortex  evolution  (via  contour  plots)  associated 
with  increasing  Reynolds  numbers.  The  corresponding  input  parameters 
are  equivalent  to  the  rectangular  channel  flow  of  Section  6.1,  and  are 
shown  in  Table  3.  A  constant  lattice/grid  resolution  of  1,000x40  is  main¬ 
tained  for  each  case. 


Figure  9.  Geometry  corresponding  to  the  flow  over  a  back-step. 


Tabie  3.  Simuiation  parameters  used  for  the  back-step  channei  flow  probiem. 


Case 

\n\etWe\.  {Uin), 

Nav-Stokes 

[m/s] 

\n\etWe\.  {Uiat) 

LB  [l/i] 

Kin.  Vise,  (v), 

Nav-Stokes 

[m2/s] 

Kin.  Vise,  (v),  LB 
[lyi] 

Re  =  200 

0.01 

0.05 

1.0  E-6 

0.01 

Re  =  400 

0.02 

0.20 

1.0  E-6 

0.02 

Re  =  800 

0.04 

0.10 

1.0  E-6 

0.005 

Re  =  1,000 

— 

0.10 

— 

0.004 
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Figure  10.  Velocity  contours  for  Re  =  200  (a),  Re  =  400  (b),  Re  =  800  (c),  and 
/?^1000  (d)  were  conducted  using  the  LBM.  The  contours  show  a  large  vortex  aft  of 
the  back-step  (as  well  as  a  secondary  circulation  zone  for  Re  =  1,000)  that  increases 
in  size  with  increasing  Reynolds  number  (lattice/grid  resolution: 100x100). 


O.IL  0.2L 


(a)  Re  =  im 


0.2L 


O.IL 


(b)  Re  =  400 


O.IL  0.2L 


(c)  Re  =  800 


(d)  1000 
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Figure  11.  Velocity  profiles  comparing  the  results  of  the  LBM  and  NS  corresponding 
to  the  flow  over  a  back-step  {Re  =  200  (a),  Re  =  400  (b)  and  /?e  =  800  (c)).  As 
indicated,  the  velocities  were  sampled  at  three  downstream  locations  corresponding 
to  VL  =  0.1,  VL  =  0.2  and  X/L  =  0.5  (lattice/grid  resolution  =  1,000x40). 


(a)  Re  =  2m 


(b)  Re  =  400 
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Figure  12.  Reattachment  distance  as  a  function  of  Reynoids  number.  Resuits  using 
the  LBM  are  compared  with  the  experiments  conducted  by  Le  and  Moin  (1994). 
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6.3.2  Results 

The  velocity  contours  conducted  with  the  LBM  are  shown  in  Figure  lo  and 
correspond  to  Re  =  200,  Re  =  400,  Re  =  800  andi?e  =  1,000.  As  shown, 
the  presence  of  an  initial  recirculation  zone  just  aft  of  the  back-step  is  evi¬ 
dent  for  each  case.  A  shear  layer  divides  this  zone  from  the  upper  stream¬ 
lines  as  shown.  The  recirculation  zone  tends  to  increase  in  size  (along  the 
X-direction)  with  increasing  Reynolds  number.  Also  shown,  for  Re  =  1000, 
is  the  presence  of  a  secondary  vortex  forming  on  the  upper  wall  with  cen¬ 
ter  at  approximately  X  =  0.3L.  As  previously  stated,  this  accumulation  of 
downstream  vortices  is  a  characteristic  of  all  laminar  incompressible  back- 
step  flows  as  long  as  the  Reynolds  number  continues  to  increase. 

The  velocity  profiles  comparing  the  results  of  the  LBM  and  Navier-Stokes 
method  are  shown  in  Figure  11.  The  results  correspond  to  Re  =  200,  Re  = 
400  and  Re  =  800.  As  shown,  the  velocity  distribution  along  the  down¬ 
stream  direction  progressively  evolves  into  a  characteristic  parabolic  pro¬ 
file  as  the  flow  becomes  fully  developed.  In  the  upstream  region  (for  X  = 
o.iL),  the  presence  of  the  back-step  is  evident  as  the  velocity  profile  be¬ 
comes  asymmetric  particularly  for  Y/H  <  0.3. 

Comparisons  between  the  LBM  and  NS  method  are  good,  particularly  at 
lower  Reynolds  numbers  (Re  =  200,  and  Re  =  400).  At  Re  =  800,  some 
disparity  between  the  methods  is  observed  at  X/L  =  0.2L,  corresponding 
to  the  circulation  region  just  short  of  the  reattachment  point,  as  seen  in 
Figure  11(c). 

Additional  comparisons  with  the  experimental  work  of  Le  et  al.  (1997) 
were  made  with  respect  to  the  reattachment  point.  As  illustrated  in  Figure 
9,  this  point  is  defined  as  the  X  location,  where  the  dividing  streamline  (lo¬ 
cated  between  the  shear  layer  and  the  recirculation  zone)  reattaches  to  the 
lower  channel  wall.  As  shown  in  Figure  12,  the  present  work  agrees  well 
with  that  of  Le  et  al.  (1997).  We  note  that  Le  et  al.  used  a  Reynolds  number 
based  on  the  maximum  velocity  (e.g.,  u^^ax  =  3u/2),  and  the  appropriate 
adjustments  were  made  here. 
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6.4  Flow  over  a  stationary  circular  cylinder 

6.4.1  Problem  description  and  setup 

The  final  benchmark  problem  involves  the  laminar  incompressible  flow 
over  a  stationary  circular  cylinder.  This  particular  flow  is  one  of  the  most 
widely  studied  in  the  CFD  literature,  and  is  well  known  for  its  ability  to  il¬ 
lustrate  the  transition  from  steady  to  unsteady  flow  behavior  as  a  function 
of  Reynolds  number.  As  shown  in  Figure  13,  at  very  low  Reynolds  number 
(Re  <  <i)  the  flow  is  steady  with  no  boundary  layer  separation.  As  Re  ap¬ 
proaches  10,  the  flow  remains  steady,  but  boundary  layer  separation  be¬ 
gins  and  evolves  to  form  two  symmetric  vortices  in  the  wake  region.  For 
Reynolds  numbers  in  excess  of  90  (but  less  than  ~104),  although  it  re¬ 
mains  laminar,  the  flow  becomes  unsteady  and  asymmetric.  Von  Karman 
vortex  shedding  is  observed  and  distinguished  by  regular  shear  layer  sepa¬ 
ration  from  the  upper  and  lower  surfaces  of  the  cylinder.  Beyond  Re  = 

-104  the  flow  becomes  turbulent. 

The  geometry  for  this  problem  is  shown  in  Figure  14.  As  shown,  a  cylinder 
of  radius  R  (R=o.im)  is  placed  6.6R  from  the  inlet  and  3R  from  the  top 
and  bottom  of  a  channel  of  length  40R.  For  boundary  conditions,  a  uni¬ 
form  inlet  velocity  is  imposed  at  the  inlet,  and  a  standard  outflow  bounda¬ 
ry  at  the  exit.  The  upper  and  lower  channel  walls  are  designated  as  no-slip 
(zero  velocity)  conditions  (implemented  with  standard  bounce-back  condi¬ 
tions  for  the  LBM  simulation). 

While  most  LBM  simulations  commonly  apply  the  standard  bounce-back 
conditions  to  the  wall  of  the  circular  cylinder,  this  action  results  (for  a  reg¬ 
ular  hexagonal  lattice)  in  an  irregular  “stair-step”  geometry.  The  effect  on 
the  flow  field  response  is  clearly  non-ideal,  as  this  irregular  surface  would 
(like  the  effect  of  dimples  on  a  golf  ball)  delay  boundary  layer  separation 
and  result  in  atypical  results.  For  this  reason,  the  IMB  method  is  invoked 
(Owen  et  al.  2011;  as  discussed  in  Chapter  6).  In  particular,  the  quantity 
(s)  that  represents  the  area  fraction  (volume  fraction  in  3D)  of  the  circular 
cylinder  with  an  intersecting  lattice  cell  is  computed  and  then  used  to 
modify  the  collision  operator. 
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Figure  13.  Illustration  of  velocity  profiles  as  a  function  of  Reynolds  number 
(http;//en.wi  kipedia.org/wiki/Reynolds_number). 


Figure  14.  Geometry  corresponding  to  the  flow  over  a  stationary  circular  cylinder. 
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For  this  work,  three  different  Reynolds  numbers  were  investigated:  ReD=i, 
ReD=io,  and  ReD=i20,  where  D  corresponds  to  the  cylinder  diameter.  The 
corresponding  input  parameters  for  both  the  NS  method  and  LBM  are 
shown  in  Table  4.  A  constant  lattice/grid  resolution  of  2000x400  is  main¬ 
tained  throughout. 

Lastly,  for  comparison  purposes  with  our  experiment,  the  Coefficient  of 
Drag  (Cd)*  is  computed  in  Eq.  66: 


*  “Coefficient  of  Drag,"  http://www.aerospaceweb.org/auestion/aerodvnamics/a0231.shtmi.  accessed 
24  Aprii  2014. 
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Cd 


Px 

XjlpU'^A 


(66) 


where  is  the  x  component  of  the  total  fluid  force  acting  on  the  cylinder, 
and  A  is  the  projected  area.  For  the  LBM,  the  total  force  (F)  over  the  cylin¬ 
der  is  computed  via  the  momentum  exchange  method,  as  shown  in  Eq.  67 
(Owen  et  al.  2011): 


F  =  #5:„s„(Lnf?i)  (67) 

where  B^,  Clf,  Sf-,  and  5^  have  been  previously  defined  in  Chapter  6. 


Table  4.  Simulation  parameters  used  for  the  back-step  channel  flow  problem. 


Case 

Inlet  Vel. 

Nav-Stokes  [m/s] 

Inlet  Vel. 

LB  [l/t] 

Kin.  Vise,  (v), 

Nav-Stokes 

[mys] 

Kin.  Vise,  (v),  LB 
p/t] 

II 

CD 

q: 

0.000005 

0.00167 

l.OE-6 

0.167 

Re=10 

0.00005 

0.0167 

1.0  E-6 

0.167 

Re=120 

0.0006 

0.2 

l.OE-6 

0.167 

6.4.2  Results 

The  velocity  contours  of  both  the  LBM  and  NS  method  are  shown  in 
Figure  15  and  correspond  to  Re  =  1,  Re  =  10,  and  Re  =  120.  As  shown,  both 
sets  of  simulation  results  are  qualitatively  equivalent,  reflecting  the 
changes  in  velocity  contours  with  respect  to  Reynolds  number  in 
accordance  with  the  expected  profiles  shown  in  Figure  13. 

Quantitative  comparisons  of  velocity  between  the  LBM  and  NS  method  are 
shown  in  Figure  16  and  correspond  to  Re  =  1  and  Re  =  10.  As  indicated,  in 
both  cases  the  velocity  distribution  at  the  downstream  location  (X/L  =  0.5) 
progressively  evolves  into  a  characteristic  parabolic  profile  as  the  flow  be¬ 
comes  fully  developed.  At  an  upstream  location  forward  of  the  cylinder 
(X/L  =  0.1),  the  presence  of  the  stationary  cylinder  is  clearly  observed  as 
the  velocity  profiles  along  Y/L  =  0.5  gradually  decrease  towards  the  point 
of  stagnation.  The  presence  of  the  expanding  wake  region  for  Re=io  is  also 
evident  at  X/L  =  0.2,  as  indicated  by  the  significant  decrease  in  velocity  at 
Y/H  =  0.5. 
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For  the  case  of  Re=i20,  the  average  drag  coefficient  (Q)  was  computed  in 
accordance  with  Eq.  66  and  67.  As  shown  in  Figure  17,  a  constant  shedding 
frequency  was  observed  after  approximately  25,000  iterations  {Uot/D  - 
12.5).  The  computed  drag  coefficient  of  0.97  compared  favorably  with  the 

experimental  results  (~0.98  - 0.99)  found  in  the  experimental  literature 

(Aerospace.org). 

Figure  15.  Velocity  contours  for  Re  =  1,  Re  =  10,  and  Re =120  show  results  from 
both  the  LBM  and  NS  method.  The  contours  reveal  the  characteristic  steady  to 
unsteady  flow  behaviors  associated  with  increasing  Reynolds  number,  and  show  the 
presence  of  Von  Karman  vortex  shedding  for  Re=  120.  (Lattice/grid  resolution: 

2000x400.) 


NS;  Re  =1 


NS;i?e=io 
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(e)  LB:7?e=  120 
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Figure  16.  Velocity  profiles  comparing  the  results  of  the  LB  and  NS  methods 
corresponding  to  the  flow  over  a  stationary  circular  cylinder  {Re  =  1,  and  /?e  =  10).  As 
indicated,  the  velocities  were  sampled  at  three  downstream  locations  corresponding 
to  X/L  =  0.1,  )Q'L  =  0.2  and  X/L  =  0.5.  (Note:  lattice/grid  resolution  =  2,000x400). 


—  X  =  0.5L  (Nav.  Stokes) 
--  X  =  0.2L  (Nav.  Stokes) 


(a)  Re=\Q 


0-0  X  =  0.1L(LBM) 

—  X  =  0.5L  (Nav.  Stokes) 


(b)  Re  =  1 
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Figure  17.  Time  history  of  the  drag  coefficient  computed  using  Eq.  66.  The  resuits 
shown  are  from  the  LBM  for  /?^120. 
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7  Summary  and  Conclusions 

This  report  documented  a  comparison/ validation  effort  accompanying  the 
development  of  a  standard  Lattice  Boltzmann  solver  which  is  to  be  coupled 
to  an  existing  ERDC  Discrete  Element  Model.  The  primary  goal  was  to 
validate  the  Lattice  Boltzmann  model  by  comparing  it  with  various 
laminar,  incompressible  flow  cases  simulated  using  a  finite  volume-based 
Navier-Stokes  solver.  Simulations  involving  four  standard  benchmark 
studies  were  analyzed:  (i)  the  flow  through  a  rectangular  channel,  (2)  the 
flow  through  a  lid-driven  cavity,  (3)  the  flow  over  a  back-step,  and  (4)  the 
flow  over  a  stationary  circular  cylinder.  Eor  these  specific  applications  and 
the  Reynolds  numbers  simulated,  the  results  showed  excellent  agreement 
between  the  two  cases.  Limitations  in  the  LBM  at  some  elevated  Reynolds 
numbers  were  due  primarily  to  compressibility  effects.  However,  even 
these  cases  showed  only  marginal  deviation  from  the  finite  volume  Navier- 
Stokes  method. 
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Appendix  A;  Derivation  of  the  Equiiibrium 
Distribution  Function 

The  equilibrium  distribution  function  can  be  derived  from  the  Maxwell- 
Boltzmann  velocity  distribution.  As  a  function  of  energy  (E),  this  can  be 
written  as  Eq.  Ai: 


/(E)  =  i4exp(-^)  (Al) 

Where  kg  is  the  Boltzmann  constant,  A  is  a  constant,  and  T  is  the  tempera¬ 
ture. 

For  a  free  gas  (particle  energy  is  solely  kinetic),  we  can  write  the  function 
in  terms  of  the  velocity  of  a  particle  as  Eq.  A2: 

f(v)  =  i4exp(-^^^)  (A2) 


Using  the  normalization  condition  (i.e.,  J  f(v)d^v  —  1),  we  obtain  the  value  for 
A,  and  f(v)  becomes  Eq.  A3: 

=  (A3, 

Rewriting  the  velocity  in  terms  of  its  mean  (u)  and  deviation  (q)  as: 

V  =  Ci-u,  and  using  the  isothermal  ideal  gas  relation:  c|  =  kgT /m,  we 
have  Eq.  A4: 


„  r  c?  1  v?-2u-Ci 

f  oc  exp - ^  ^^P - ^ — ' 


(A4) 


Using  the  Taylor  series  expansion  (e^  =  l  +  x  +  —  +■■•)  and  dropping 
terms  of  order  three  and  above,  gives  Eq.  A5: 


exp 


v?-2u-c{ 
2c?  . 


U-Cj  ,  u 


1  +  ^  + 


2c? 


2c 


(A5) 
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Since  along  with  the  collision  operator  (H)  the  equilibrium  function  must 
preserve  both  mass  and  momentum,  as  shown  in  Eq.  A6  and  Eq.  A7: 


P  =  Si =  lift 

(A6) 

pu  =  =  HiCifi 

(A7) 

It  can  be  shown  that  the  constant  of  proportionality  (K)  is  simply  p 
(Viggen  2009).  This  gives  us  the  final  generalize  form  of  the  equilibrium 
function  as  Eq.  A8: 


U-Cl  (u-Cj)^ 

cl  2cl 
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